function [ll, est] = measurement(xparam, hhparam, psi0, PP0, y, y5, inflation, forecast)

[T, N] = size(y); Ny = 1+2*N; Nxi = 3+5*N;

v_eta = exp(xparam(1)); 
v_nu1 = exp(xparam(2)); 
v_nu2 = exp(xparam(3)); 
v_e   = exp(xparam(4)); 

hh.v_w = exp(hhparam(1:N));     
hh.v_e = exp(hhparam(N+1:2*N)); 
hh.theta = hhparam(2*N+1:3*N);  hh.theta(1)=0;

hh.load = zeros(N,1);

if forecast == 0

    psi0(2:4:2+4*N) = hhparam(3*N+1:4*N+1); 
    psi0(3:4:3+4*N) = hhparam(4*N+2:5*N+2);

end

H = zeros(Ny, Nxi); H(1,1) = 1;

R = diag([v_e; hh.v_e; hh.v_e]);

Q = zeros(Nxi); Q(1,1) = v_eta; Q(2,2) = v_nu1; Q(3,3) = v_nu2;

F = zeros(Nxi); F(2:3,2:3) = [1,0;0,1];

hh.K = NaN(3,1,T,N); hh.P = NaN(3,3,T+1,N);

for n = 1:N
    
    H(1+n, 3+4*(n-1)+2) = 1; H(1+N+n, 3+4*N+n) = 1;
    
    Q(4*n,1) = v_eta; Q(1,4*n) = v_eta;
    Q(3+4*(n-1)+1,3+4*(n-1)+1) = v_eta + hh.v_w(n) + hh.load(n)^2;  
    
    hh.P(:,:,1,n) = diag([1;1;1]);  % HH prior
    
    for ii = n+1:N
        Q(3+4*(n-1)+1,3+4*(ii-1)+1) = hh.load(n)*hh.load(ii);
        Q(3+4*(ii-1)+1,3+4*(n-1)+1) = hh.load(n)*hh.load(ii);
    end
end

PP = NaN(Nxi, Nxi,T+1); PP2 = NaN(Nxi, Nxi,T+1); KK = NaN(Nxi, Ny,T); FF = NaN(Nxi, Nxi,T);
psi_1 = NaN(Nxi,1,T+1); psi_2 = NaN(Nxi,1,T); u = NaN(Ny,1,T); 
cov_u = NaN(Ny, Ny, T); 

PP(:,:,1) = PP0;
psi_1(:,1,1) = psi0; % psi(0|0)

ll = 0;
cc = 0.1;

for t = 1:T

    Rt = R; Ht = H;

    f_phi = 1/(1+cc*exp(-psi_1(3,1,t)));
    df_phi = cc*exp(psi_1(3,1,t))/(cc+exp(psi_1(3,1,t)))^2;
    
    FF(:,:,t) = F;
    FF(1,1:3,t) = [f_phi, 1, psi_1(1,1,t) * df_phi];
    
    psi_2(:,1,t) = psi_1(:,1,t); 
    psi_2(1,1,t) = psi_1(2,1,t) + psi_1(1,1,t) * f_phi;  % psi(t|t-1)
    
    for n = 1:N
        s_it = psi_1(3+4*(n-1)+1,1,t);
        pi = psi_1(3+4*(n-1)+2,1,t);
        alpha = psi_1(3+4*(n-1)+3,1,t); 
        phi = 1/(1+cc*exp(-psi_1(3+4*n,1,t)));
        d_phi = cc*exp(psi_1(3+4*n,1,t))/(cc+exp(psi_1(3+4*n,1,t)))^2;
        
        [hh.K(:,:,t,n) , hh.P(:,:,t+1,n)] = hh_K(v_eta, v_nu1, v_nu2, hh.v_w(n), hh.P(:,:,t,n), phi, pi, d_phi);
        K_a = hh.K(2,1,t,n);
        K_phi = hh.K(3,1,t,n);
        K_pi = hh.K(1,1,t,n);

        psi_2(3+4*(n-1)+1,1,t) = psi_1(2,1,t) + psi_1(1,1,t) * f_phi + hh.theta(n); % s(t|t-1)
        psi_2(3+4*(n-1)+2,1,t) = alpha + phi * pi; % SR + hh.theta(n)
        psi_2(3+4*N+n,1,t) = (alpha)/(1-phi); % LR
        
        FF(3+4*(n-1)+1,1:3,t) = [f_phi, 1, psi_1(1,1,t) * df_phi];
        FF(3+4*(n-1)+2:3+4*n,3+4*(n-1)+1:3+4*n,t) = ...
            [ K_a + K_phi*pi*d_phi + K_pi*phi + 2*K_phi*K_pi*(s_it-pi)*d_phi,...
             -K_a + phi - 2*K_phi*pi*d_phi + K_phi*s_it*d_phi - K_pi*phi - 2*K_phi*K_pi*(s_it-pi)*d_phi, ...
             1, ...
             pi + K_pi*(s_it - pi);
             K_a,   -K_a,   1,  0;
             K_phi, -K_phi, 0,  1];
        FF(3+4*N+n, 3+4*n-1:3+4*n,t) = [1/(1-phi), d_phi*alpha/(1-phi)^2];

        % handle missing values in LR expectations
        if isnan(y5(t,n))
            
            y5(t,n) = 0;
            Rt(1+N+n, 1+N+n) = 1;
            Ht(1+N+n, 3+4*N+n) = 0;       

        end

        if isnan(y(t,n))
            
            y(t,n) = 0;
            Rt(1+n, 1+n) = 1;
            Ht(1+n, 3+4*(n-1)+1) = 0;       

        end

    end
    
    P_t = FF(:,:,t) * PP(:,:,t) * FF(:,:,t)' + Q; PP2(:,:,t) = P_t;                  % P_t = P{t|t-1}
    cov_u(:,:,t) = Ht * P_t * Ht' + Rt;                   %inv_v_u = inv(cov_u(:,:,t));
    KK(:,:,t) = P_t * Ht' / cov_u(:,:,t);     
    PP(:,:,t+1) = P_t - P_t * Ht' / cov_u(:,:,t) * Ht * P_t;    % PP(t+1) = P{t|t} 
    u(:,1,t) = [inflation(t),y(t,:),y5(t,:)]' - Ht*psi_2(:,1,t);
    psi_1(:,1,t+1) = psi_2(:,1,t) + KK(:,:,t) * u(:,1,t); %psi_1(t+1) = psi(t|t)
    
    %if t>1
    ll = ll + 0.5 * log(det(cov_u(:,:,t))) + 0.5 * u(:,1,t)' / cov_u(:,:,t) * u(:,1,t);
    %end
    
end

% shat = xi{t|t}, sfor = xi{t|t-1}
% sig = P{t|t}, omega = P{t|t-1}

est.sig = PP;
est.omega = PP2;
est.shat = psi_1;
est.sfor = psi_2;
est.pi_re = psi_2(1,1,:);

est.cov = cov_u;
est.hh_P = hh.P;
est.hh_K = hh.K;
est.KK = KK;
est.u = reshape(u,Ny,T);
est.FF = FF;
est.H = H;


end